Motivation - Deutschlands größte Klimasünder

• Spiegel Artikel am 16.11.2017 aus Anlass der Jamaika Gespräche

Gliederung

• Quellen für Polygonzüge (Staaten, Gemeinden, PLZ Bereiche etc.)
• Quellen für inhaltliche Daten
• Pakete zur Erstellung thematischer Karten (bspw. Pakete maptools, sp, tmap)
• Verknüpfung von Daten
• Beispiele für die Darstellung in Karten

Quellen für Polygonzüge

Hello World

library(maps)
map()

Das Paket maps - etwas detailierter

Grenzen sind recht grob:

map("world", "Germany")

Das Paket maps - Mehr Information

• Nur für manche Staaten bekommt man Umkreise für Einheiten unterhalb der Staatsgrenze (bspw. Frankreich, USA).

data(world.cities)
map("france")
map.cities(world.cities,col="blue")

Das Paket maptools

• Das Paket maptools hat intuitivere Bedienung, zudem können Shapefiles verarbeitet werden.

library(maptools)
data(wrld_simpl)
plot(wrld_simpl,col="royalblue")

Was sind shapefiles (.shp)?

• Das Dateiformat Shapefile ist ein ursprünglich für die Software ArcView der Firma ESRI entwickeltes Format für Geodaten. (Quelle: Wikipedia)

head(wrld_simpl@data)

length(wrld_simpl)

## [1] 246

nrow(wrld_simpl@data)

## [1] 246

Einzelne Elemente des Datensatzes plotten

ind <- which(wrld_simpl$ISO3=="DEU")

plot(wrld_simpl[ind,])

wrld_simpl@data[ind,]

##     FIPS ISO2 ISO3  UN    NAME  AREA  POP2005 REGION SUBREGION   LON   LAT
## DEU   GM   DE  DEU 276 Germany 34895 82652369    150       155 9.851 51.11

Das R-Paket choroplethrMaps

library(ggplot2)
library(choroplethrMaps)
data(country.map)
ggplot(country.map, aes(long, lat, group=group)) + geom_polygon()

Eine Karte für die USA

data(state.map)
ggplot(state.map,aes(long,lat,group=group))+geom_polygon()

Andere Quellen für Shapefiles - Das Paket raster

• Für Polygonzüge unterhalb der Staatsgrenzen ist Global Administrative Boundaries eine gute Quelle.
• Vor allem weil es auch eine API gibt, die man mit dem Paket raster nutzen kann.

library(raster)
LUX1 <- getData('GADM', country='LUX', level=1)
plot(LUX1)

Daten für das Luxemburg Beispiel

head(LUX1@data)

Shapefiles bei Eurostat

BKG - Quelle für Kreise in Deutschland

• Umrisse von 402 Kreisen in Deutschland
• Quelle: Bundesamt für Kartographie und Geodäsie (BKG)
• Karten gibt es auch für Bundesländer und Gemeinden

library(maptools)
krs <- readShapePoly("vg250_ebenen/vg250_krs.shp")
plot(krs)

head(krs@data$RS)

## [1] 03401 03458 09473 05962 10046 05916
## 402 Levels: 01001 01002 01003 01004 01051 01053 01054 01055 01056 ... 16077

Die Kreise für Baden-Württemberg

• Systematik hinter dem Amtlichen Gemeindeschlüssel (AGS) bzw. Regionalschlüssel (RS) bei Wikipedia

BLA <- substr(krs@data$RS,1,2)
plot(krs[BLA=="08",])

Shapefiles für Wahlkreise

Ortsnetzbereiche

Quelle: Bundesnetzagentur

onb <- readShapePoly("onb_grenzen.shp")

head(onb@data)

Karte der Vorwahlbereiche

Einen größeren Vorwahlbereich ausschneiden

vwb <- as.character(onb@data$VORWAHL)
vwb1 <- substr(vwb, 1,2)
vwb7 <- onb[vwb1=="07",]
plot(vwb7)

vwb <- as.character(onb@data$ONB_NUMMER)
vwb1 <- substr(vwb, 1,1)
vwb7 <- onb[vwb1=="7",]
plot(vwb7)

Das Paket rgdal

• Postleitzahlenbereiche - http://arnulf.us/PLZ

library(rgdal)

## OGR data source with driver: ESRI Shapefile 
## Source: "post_pl.shp", layer: "post_pl"
## with 8270 features
## It has 3 fields

library(rgdal)
PLZ <- readOGR ("post_pl.shp","post_pl")

PLZ-Bereiche in Stuttgart

SG <- PLZ[PLZ@data$PLZORT99=="Stuttgart",]
plot(SG,col="chocolate1")

PLZ-Bereiche in Berlin

BE <- PLZ[PLZ@data$PLZORT99%in%c("Berlin-West","Berlin (östl. Stadtbezirke)"),]
plot(BE,col="chocolate2",border="lightgray")

Zwischenfazit - Quellen für Polygonzüge

• In einigen R-Paketen sind Polygonzüge hinterlegt
• Umgang mit den Daten unterscheidet sich
• Am sinnvollsten ist das arbeiten mit Shapefiles (.shp)
• Neben den R-Paketen gibt es zahlreiche weitere Quellen für Polygonzüge.

Thematische Karten mit R erstellen

Thematische Karten erzeugen - das Paket sp

library(sp)
spplot(wrld_simpl,"POP2005")

Andere Einfärbungen wählen - Das Paket colorRamps

library(colorRamps)
spplot(wrld_simpl,"POP2005",col.regions=blue2red(100))

Es gibt auch noch mehr Farbverläufe im Paket colorRamps

• blue2green, blue2yellow

spplot(wrld_simpl,"POP2005",col.regions=matlab.like(100))

Thematische Karten mit dem Paket choroplethr

library(choroplethr)
data(df_pop_state)

head(df_pop_state)

Eine thematische Karte mit choroplethr erstellen

state_choropleth(df_pop_state)

Nur drei Staaten darstellen

state_choropleth(df_pop_state,
                 title      = "2012 Population Estimates",
                 legend     = "Population",
                 num_colors = 1,
                 zoom       = c("california", "washington", "oregon"))

Eine Karte der US Counties

data(df_pop_county)
county_choropleth(df_pop_county)

Choroplethen Länder

data(df_pop_country)
country_choropleth(df_pop_country,
              title      = "2012 Population Estimates",
              legend     = "Population",
              num_colors = 1,
              zoom       = c("austria","germany",
                             "poland", "switzerland"))

Weltbank Daten

library(WDI) 
WDI_dat <- WDI(country="all", indicator=c("AG.AGR.TRAC.NO","TM.TAX.TCOM.BC.ZS"),
    start=1990, end=2000)

• Es gibt auch eine Funktion WDIsearch mit der man nach Indikatoren suchen kann

head(WDI_dat)

Weltkarte mit den Weltbank Daten

choroplethr_wdi(code="SP.DYN.LE00.IN", year=2012,
                title="2012 Life Expectancy")

Eurostat Daten

Sie können eine Statistik der Sparquote bei Eurostat downloaden.

http://ec.europa.eu/eurostat/web/euro-indicators/peeis

library(xlsx)
HHsr <- read.xlsx2("data/HHsavingRate.xls",1)

Zensus Ergebnisse

Zensus Atlas

https://ergebnisse.zensus2011.de/

Zensus Datenbank

Zensus Gemeindeergebnisse

zen <- read.csv2("data/Zensus_extract.csv")
# Personen mit eigener Migrationserfahrung
# mit beidseitigem Migrationshintergrund
zen2 <- data.frame(Personen_Mig=zen[,which(zen[9,]==128)],
                   Personen_Mig_bs=zen[,which(zen[9,]==133)])

datahub.io

Weltkulturerbestätten

url <- "https://raw.githubusercontent.com/Japhilko/
GeoData/master/2015/data/whcSites.csv"

whcSites <- read.csv(url) 

Exkurs - OpenStreetMap Projekt

OpenStreetMap.org ist ein im Jahre 2004 gegründetes internationales Projekt mit dem Ziel, eine freie Weltkarte zu erschaffen. Dafür sammeln wir weltweit Daten über Straßen, Eisenbahnen, Flüsse, Wälder, Häuser und vieles mehr.

http://www.openstreetmap.de/

OSM Mapper

Export von OpenStreetMap Daten

<www.openstreetmap.org/export>

OpenStreetMap - Map Features

Overpass Turbo

https://overpass-turbo.eu/

Query Overpass

• Overpass hat eine eigene Abfragesprache

node
  [amenity=bar]
  ({{bbox}});
out;

Zwischenfazit - Quellen für inhaltliche Daten

• Gerade auf Staatenebene gibt es sehr viele Daten
• Oftmals hat man aber eigene Daten, die man in einer Karte darstellen möchte
• Frage ist wie sich die inhaltlichen Daten mit den Polygonzügen verbinden lassen

Verknüpfung von Daten

Quelle: Geographic Information Systems and Remote Sensing

Daten verbinden - Beispiel Eurostat Daten

• Wir hatten vorhin schon die Haushaltssparrate geladen
• Dabei handelt es sich um einen von den Ausgewählte Wichtigsten Europäische Wirtschaftsindikatoren bei Eurostat
• Der Datensatz wrld_simpl aus dem Paket maptools wird verwendet
• Die Daten werden mit Hilfe des Ländernamens miteinander verknüpft

ind <- match(HHsr$geo,wrld_simpl@data$NAME)
ind <- ind[-which(is.na(ind))]

• Der wrld_simpl Datensatz wird auf Europa eingeschränkt

EUR <- wrld_simpl[ind,]

• Informationen zur Haushaltssparrate werden hinzugefügt

EUR@data$HHSR_2012Q3 <- as.numeric(as.character(HHsr[-(1:2),2]))
EUR@data$HHSR_2015Q2 <- as.numeric(as.character(HHsr[-(1:2),13]))

Karte mit Eurostat Indikator Household Saving Rate

• Mit dem Befehl spplot aus dem Paket sp ist es möglich zwei Karten nebeneinander darzustellen:

spplot(EUR,c("HHSR_2012Q3","HHSR_2015Q2"))

Daten verbinden - Beispiel Bäckereien in Berlin

• Quelle für die folgenden Daten ist:

OSM als Datenquelle

• Zum Download habe ich die Overpass API verwendet

(load("data/info_bar_Berlin.RData"))

## [1] "info"

Verwendung des Pakets gosmd

devtools::install_github("Japhilko/gosmd")

library("gosmd")
pg_MA <- get_osm_nodes(object="leisure=playground","Mannheim")
pg_MA <- extract_osm_nodes(pg_MA,value='playground')

Matching

tab_plz <- table(info_be$addr.postcode)

ind <- match(BE@data$PLZ99_N,names(tab_plz))
ind

##   [1]  1  2  3  4  5  6  7  8 NA  9 NA NA NA NA NA 10 11 12 NA 13 14 15 16
##  [24] 17 18 19 20 21 22 23 24 25 NA 26 27 28 29 NA NA NA NA 30 NA 31 32 33
##  [47] 34 35 NA NA 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 NA 52 53
##  [70] NA 54 55 NA NA NA 56 57 58 59 60 NA NA NA NA NA 61 NA NA NA 62 NA NA
##  [93] NA NA NA NA NA NA NA 63 NA NA 64 NA 65 NA NA NA 66 NA NA NA NA 67 NA
## [116] NA 68 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [139] NA 69 70 NA 71 72 73 74 75 NA 76 NA NA NA NA NA NA NA NA NA NA NA NA
## [162] 77 NA 78 79 NA NA NA NA 80 NA NA NA NA 81 NA 82 83 84 NA NA NA NA NA
## [185] NA NA NA 85 NA NA

Daten anspielen

BE@data$num_plz <- tab_plz[ind]

Das Paket tmap

library(tmap)

• mit qtm kann man schnell eine thematische Karte erzeugen

BE@data$num_plz[is.na(BE@data$num_plz)] <- 0
qtm(BE,fill = "num_plz")

Mehr Informationen einbinden

• Der folgende Datensatz ist eine Kombination aus den vorgestellten PLZ-Shapefiles und OSM-Daten die über Overpass heruntergeladen wurden:

load("data/osmsa_PLZ_14.RData")

OSM-Daten - Bäckereien in Stuttgart

qtm(PLZ_SG,fill="bakery")

In welchem PLZ Bereich sind die meisten Bäckereien

kable(PLZ_SG@data[which.max(PLZ_SG$bakery),c("PLZ99","lat","lon","bakery")])

Das R-Paket RDSTK

• Data Science Toolkit API

library("RDSTK")

Die Daten für Stuttgart

PLZ_SG <- PLZ[PLZ@data$PLZORT99=="Stuttgart",]

Eine Karte der Flächenbedeckung erstellen

• Daten von European Commission Land Resource Management Unit Global Land Cover 2000.

qtm(PLZ_SG,fill="land_cover.value")

Die Höhe in Stuttgart

• Daten von NASA and the CGIAR Consortium for Spatial Information .

qtm(PLZ_SG,fill="elevation.value")

Graphiken Stadtleben Stuttgart - das Paket ggmap

devtools::install_github("dkahle/ggmap")
install.packages("ggmap")

Eine erste Karte mit ggmap erzeugen

library(ggmap)

qmap("Stuttgart")

Karte für einen ganzen Staat

qmap("Germany")

• Wir brauchen ein anderes zoom level

Ein anderes zoom level

• level 3 - Kontinent
• level 10 - Stadt
• level 21 - Gebäude

qmap("Germany", zoom = 6)

Karte für eine Sehenswürdigkeit

• ggmap - maptype satellite

WIL <- qmap("Wilhelma",zoom=20, maptype="satellite")
WIL

ggmap - maptype satellite zoom 20

qmap('Stuttgart Hauptbahnhof', zoom = 15, maptype="hybrid")

Eine terrain Karte

qmap('Stuttgart Fernsehturm', zoom = 14,
 maptype="terrain")

ggmap - maptype watercolor

qmap('Stuttgart', zoom = 14,
 maptype="watercolor",source="stamen")

ggmap - source stamen

qmap('Stuttgart', zoom = 14,
 maptype="toner",source="stamen")

ggmap - maptype toner-lite

qmap('Stuttgart', zoom = 14,
 maptype="toner-lite",source="stamen")

# ggmap - maptype toner-hybrid

qmap('Stuttgart', zoom = 14,
 maptype="toner-hybrid",source="stamen")

ggmap - maptype terrain-lines

qmap('Stuttgart', zoom = 14,
 maptype="terrain-lines",source="stamen")

Geokodierung

Geocoding (…) uses a description of a location, most typically a postal address or place name, to find geographic coordinates from spatial reference data …

Wikipedia - Geocoding

library(ggmap)
geocode("Stuttgart")

##        lon      lat
## 1 9.182932 48.77585

Reverse Geokodierung

Reverse geocoding is the process of back (reverse) coding of a point location (latitude, longitude) to a readable address or place name. This permits the identification of nearby street addresses, places, and/or areal subdivisions such as neighbourhoods, county, state, or country.

Quelle: Wikipedia

revgeocode(c(48,8))

## [1] "Unnamed Road, Somalia"

Die Distanz zwischen zwei Punkten

mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart")

##                    from                     to    m    km   miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 3136 3.136 1.94871     488
##    minutes     hours
## 1 8.133333 0.1355556

mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart",mode="walking")

##                    from                     to    m    km    miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 2505 2.505 1.556607    1874
##    minutes     hours
## 1 31.23333 0.5205556

Eine andere Distanz bekommen

mapdist("Marienplatz Stuttgart","Hauptbahnhof Stuttgart",mode="bicycling")

##                    from                     to    m    km    miles seconds
## 1 Marienplatz Stuttgart Hauptbahnhof Stuttgart 2722 2.722 1.691451     486
##   minutes hours
## 1     8.1 0.135

Geokodierung - verschiedene Punkte von Interesse

POI1 <- geocode("B2, 1 Mannheim",source="google")
POI2 <- geocode("Hbf Mannheim",source="google")
POI3 <- geocode("Mannheim, Friedrichsplatz",source="google")
ListPOI <-rbind(POI1,POI2,POI3)
POI1;POI2;POI3

##        lon      lat
## 1 8.462844 49.48569

##        lon      lat
## 1 8.469879 49.47972

##        lon      lat
## 1 8.475754 49.48304

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),
data = ListPOI)

Punkte in der Karte

MA_map +
geom_point(aes(x = lon, y = lat),col="red",
data = ListPOI)

ggmap - verschiedene Farben

ListPOI$color <- c("A","B","C")
MA_map +
geom_point(aes(x = lon, y = lat,col=color),
data = ListPOI)

ggmap - größere Punkte

ListPOI$size <- c(10,20,30)
MA_map +
geom_point(aes(x = lon, y = lat,col=color,size=size),
data = ListPOI)

Eine Route von Google maps bekommen

from <- "Mannheim Hbf"
to <- "Mannheim B2 , 1"
route_df <- route(from, to, structure = "route")

Mehr Information

http://rpackages.ianhowson.com/cran/ggmap/man/route.html

Eine Karte mit dieser Information zeichnen

qmap("Mannheim Hbf", zoom = 14) +
  geom_path(
    aes(x = lon, y = lat),  colour = "red", size = 1.5,
    data = route_df, lineend = "round"
  )

Wie fügt man Punkte hinzu

• Nutzung von geom_point

• Question on stackoverflow

http://i.stack.imgur.com

Cheatsheet

• Cheatsheet zu data visualisation

https://www.rstudio.com/

Das Paket ggmap

library(ggmap)
lon_plz <- PLZ_SG@data[which.max(PLZ_SG$bakery),"lon"]
lat_plz <- PLZ_SG@data[which.max(PLZ_SG$bakery),"lat"]
mp_plz <- as.numeric(c(lon_plz,lat_plz))
qmap(location = mp_plz,zoom=15)

Das Paket osmar

library(osmar) 

src <- osmsource_api()
gc <- geocode("Stuttgart-Degerloch")
bb <- center_bbox(gc$lon, gc$lat, 800, 800)
ua <- get_osm(bb, source = src)
plot(ua)

Gebäude in diesem Ausschnitt plotten

bg_ids <- find(ua, way(tags(k=="building")))
bg_ids <- find_down(ua, way(bg_ids))
bg <- subset(ua, ids = bg_ids)
bg_poly <- as_sp(bg, "polygons")  
plot(bg_poly)

Gebäude und Straßen im Ausschnitt

Schlussfolie

Vielen Dank für die Aufmerksamkeit

Erweiterungsmöglichkeiten

• Interaktive Karten mit dem Paket leaflet

http://rpubs.com/Japhilko82/Rleaflet

Karte auf Carto zur Energieerzeugung